Three-dimensional relativistic particle-in-cell hybrid 
code based on an exponential integrator 
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Abstract — In this paper we present a new three dimensional 
(3D) full electromagnetic relativistic hybrid plasma code H-VLPL 
(hybrid virtual laser plasma laboratory). The full kinetic particle- 
in-cell (PIC) method is used to simulate low density hot plasmas 
while the hydrodynamic model applies to the high density 
cold background plasma. To simulate the linear electromagnetic 
response of the high density plasma, we use a newly developed 
form of an exponential integrator method. It allows us to simulate 
plasmas of arbitrary densities using large time steps. The model 
reproduces the plasma dispersion and gives correct spatial scales 
like the plasma skin depth even for large grid cell sizes. We 
test the hybrid model validity by applying it to some physical 
examples. 

I. Introduction 

It is well accepted now that the particle-in-Cell (PIC) codes 
provide the most detailed description of plasmas and are the 
key computational tools in the study of relativistic laser-plasma 
interactions [1], [2]. Large full 3D parallel electromagnetic 
simulation codes like VLPL [3], OSIRIS [4], VORPAL [5], 
OOPIC [6], and others contributed remarkably in our under- 
standing of the complex laser-plasma physics. Because these 
codes provide the most detailed plasma description, they are 
computationally expensive. As a result one continues to look 
for new algorithms and simulation techniques to cope with 
challenges of the laser-plasma physics. 

One of the reasons why the classical explicit PIC methods are 
computationally extremely expensive is that they have to re- 
solve the fundamental plasma frequency ujp = y^^Trn^e^Jrn^, 
which is the frequency of the plasma electrostatic oscillations. 
Therefore, they are efficient only when applied to low density 
plasmas. 

At the same time, there is a number of important applications 
where lasers interact with high density plasmas, e.g., the 
studies of electron propagation through solid density targets 
and the resulting target normal sheath acceleration (TNSA) 
[7]. The solid state density plasma densities are in the range 
of 100 - 1000 Ur, where 



Tie = muj^ /4:7Te'^ 



(1) 



is the critical plasma density. Here, m is the electron mass, 
— e is its charge, and uj is the laser frequency. For the 1 /im 
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wavelength laser the critical electron density is 10^^ 1/cc. 
Other important applications include the Fast Ignition (FI) 
physics in the Inertial Confinement Fusion (ICF) studies [8]. 
The FI plasma has a density of the 1000 times compressed 
solid hydrogen, i.e., of the order of 10^ ric. Hence, the 
applicability of the classical PIC codes in this density range 
is facing severe difficulties. In this situation, one is forced 
to look for a more efficient numerical method to challenge 
those ultra-high densities. One of the possibilities is to include 
hydrodynamic description of the high density plasma in the 
fully kinetic PIC code. 

In last years PIC-hydrodynamic hybrid techniques have 
emerged as an efficient solution to large scale ultra high- 
density plasma simulations, e.g., FI physics, solid state density 
plasma interactions, high charge, high energy ion generations 
etc [9]-[ll]. Most of these codes work in the Darwin approxi- 
mation and thus exclude the electromagnetic wave propagation 
completely. They also exclude electrostatic waves keeping 
the collisional magnetohydrodynamics (MHD) only. Further, 
an implicit electrostatic particle-fluid hybrid plasma code has 
been developed by Rambo and Denavit [12], which has been 
used to study interpenetration and ion separation in colliding 
plasmas [13]. There is also the implicit electromagnetic PIC 
code LSP [14]. This code uses an implicit global scheme 
which overcomes such restrictions of the time- step. The LSP 
code also employs a field solver based on an unconditionally 
Courant- stable algorithm [15] for electromagnetic calculations. 

Recently, we have presented a ID version of the code Hybrid 
Virtual Laser Plasma Laboratory (H-VLPL) [16] that unites 
a hydrodynamic model for overdense plasmas and the full 
kinetic description of hot low-density electrons and ions. In 
this code, the linear plasma response was simulated using an 
implicit scheme. The implementation involved the solution of 
linear systems, which have been done in a very efficient way 
using the Schur complement. 

Unfortunately, the efficiency of the implicit H-VLPL code 
drops significantly if we extend the code from ID to 3D. 
Therefore, we introduce a new 3D version of the code that 
is based on a different approach. Instead of using an implicit 
method, we employ a special variant of an exponential in- 
tegrator [17] to model the high frequency plasma response. 



Exponential integrators are methods which make use of matrix 
functions related to the matrix exponential of the Jacobian of 
the differential equation. Here we consider a modification of 
the mollified impulse method [18], which has been proposed 
for molecular dynamics simulations. 

The mollified impulse method is motivated from a splitting 
approach. Variants of splitting methods are widely used for 
problems acting on different time scales, see [19]. For our 
hybrid model it turns out that due to the high density of 
the plasma, the highest frequencies stem from a multplication 
operator, which acts only locally on each grid point. Frequen- 
cies arising from the Maxwellian part are much lower and 
can be handled explicitly as in the PIC code. This allows to 
implement the new mollified impulse method by evaluating 
matrix functions of diagonal matrices only. Obviously, this 
is much more efficient than the solution of linear systems 
resulting from a 3D discretization. 

To illustrate the performance of the new method, we apply it to 
a few test physics examples. We check the correct dispersion of 
electromagnetic waves in the hybrid plasma and also compare 
the numerical skin length with the analytical expressions. 
As a more complicated test, we apply the new 3D code to 
the laser- solid interaction and acceleration of protons in the 
TNSA regime. Additionally, we perform a feasibility study of 
simulations of the Weibel instability, which occurs in the FI 
scenario. 

The paper is organized as follows. First, we describe the 
full hybrid method in Section II. Then, we briefly explain 
the numerical scheme in Section III and finally, we test the 
new code H-VLPL on some well-known physical examples in 
Section IV. 



II. Hybrid model 



In this section we define the physical model we use to simulate 
the plasma. We use the full kinetic description with the usual 
PIC macroparticles for low density hot electrons and ions. 
The high density cold background plasma is then described 
hydrodynamically. Since the spatial locations of the kinetic 
particles and the hydrodynamic parts may overlap, we have to 
add currents generated by all species in the same grid cells. 

The equations for the fields and particle momenta read 



— =cVxB-47r^j„ i = , 

dt ^ 

at c 



i, h 



(2a) 

(2b) 
(2c) 
(2d) 



where 



j^ = qini\i, p^ = m^7^v^, 



l£ = Wl 



{rriicY 



The index £ = e^i^h denotes electrons, ions, and hybrid par- 
ticles, respectively. E and B denote the electric and magnetic 
field vectors, j denotes the current density, p is the momentum 
and n the number density of particles. 

In the momentum equation (2c) we have neglected the non- 
linear part of the Lorentz force v x B/c, because we assume 
that the velocity of the electron part in the cold background 
plasma is small, v <C c. This assumption, however, limits the 
cold plasma response to the linear one. 

III. Numerical algorithm 

For simplicity, we rewrite the equations in dimensionless 
variables, t = uj^t and x = k^x, where ujq denotes the laser 
frequency and ko = ujq/c. The new variables are then 
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where i = e,i. 













In the following, we omit the tildes, neglect the Lorentz force 
and consider hybrid particles only. We can then write p = P/^. 
Using these simplification, Eq. (2) reads 







f = VxB + .?p 


(3a) 






5B „ ^ 

dt= ^^^ 


(3b) 






dv 


(3c) 


where uj^ = 


rih 
Ih ' 







The problem is considered in three space dimensions. We 
solve the equations on a staggered grid and approximate the 
spatial derivatives with centered finite differences using the 
Yee scheme [20]. 

For the time discretization we will use the following splitting 
of the vector fields 
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(4) 



If ujp is constant over a time step, the exact solution of 
each of the three differential equations can be computed very 



efficiently, in particular without solving any linear system. A 
symmetric splitting yields the following scheme 
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Although this splitting method is of classical order two (since 
it is a symmetric scheme), it suffers from resonances, which 
arise in E, B, and p if the density becomes large. In fact, 
the errors of this scheme are of order zero for certain time 
steps. We illustrate this effect by simulating a ID plane wave. 
The incoming laser pulse is modeled via inhomogeneous, time 
dependent Dirichlet boundary conditions and zero as initial 
values. A spatial grid size of 0.5 for x G [0, 200] is used. 
The hybrid density is set to rih = lO^ric and the system is 
integrated over the time interval [0,200]. The blue curves in 
Fig. 1 shows the errors in Ey, B^ and Py of the standard 
splitting method (5) as a function of the time step size, while 
the red line corresponds to a second-order error behavior. To 
improve the presentation, we only show the interval [0.25, 0.5] 
for the time steps, but we would like to emphasize that the 
same effects have been obtained for much smaller time steps 
as well. 
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Fig. 1. Error in Ey, Bz and py plotted over the step size (blue) for a 
straightforward integration of Eq. (4). The red line shows the expected order 
two. 



Similar resonance effects have also been observed for multiple 
time stepping schemes in molecular dynamics simulations 
[21] and for numerical methods for solving second-order 
differential equations [18], [22], [23]. Motivated by these 
papers, we suggest to apply filter functions and averaging 



operators to the Maxwellian part and modify the standard 
splitting method (5) in the following way 
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Fig. 2 shows the same numerical test as for Fig. 1 with 
filter functions (j){x) = ip{x) = sinc(a;/2), where sinc(a;) := 
sin(a;)/a; and time steps t e [0.1,0.5]. The resonances 
have been eliminated completely and order two is achieved 
for arbitrarily large densities. The theoretical properties of 
the numerical method including a detailed error analysis is 
currently investigated and will be reported elsewhere. 
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Fig. 2. Error in Ey, Bz and py plotted over the step size (blue) with 
(f)(^x) = ip{x) = sine |. The red line indicates the expected order two. 



IV. Benchmark with physical processes 

The numerical integrator which was described in the previous 
sections has been implemented into the VLPL code as a three- 
dimensional, parallelized version, and is now operational. 
In order to examine its accuracy and reliability, we have 
benchmarked it with a variety of physical processes. 
First, we check if it correctly models laser propagation 
through linearly dispersive plasma as well as reflection from 
overdense plasma. Second, we verify the conservation of the 
total energy of the system by the hybrid algorithm. Third, 
our code is applied to the Target Normal Sheath Acceleration 
(TNSA) process, which would have been very difficult to treat 
just using PIC means since it uses materials of solid state 
density. We check if our hybrid integrator correctly describes 



the exponential decay of a wave in overdense plasma. Finally, 
we show its applicability to study the Weibel instability. 



A. Reflection of an incident pulse 

As the simplest test one can imagine, we will show that our 
integrator accurately models refraction in underdense plasma 
and reflection from overdense plasma. First, we set up a 
plasma slab of 0.85nc density (1) and send a 26/s Gaussian 
laser pulse through it. As the pulse hits the surface of the 
purely hybrid plasma, a part of the wave is transmitted while 
a significant reflection also occurs. 
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Fig. 3. Snapshots of the simulation setup taken with an interactive VR 
visualization software, which is a part of H-VLPL and currently under 
development. The left picture shows the laser pulse (isosurfaces of fixed 
positive and negative electric field amplitudes) as it enters the hybrid plasma. 
The right picture demonstrates the dispersive effect. 
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Fig. 4. Plot of the total energy, the energy of the electromagnetic field and 
of the hybrid plasma versus time. 



During the laser propagation in vacuum, the energy stays 
constant except for small fluctuations within the order of 
magnitude of the machine precision. When the pulse hits the 
overdense hybrid plasma surface, it is reflected, as can be seen 
at the time of 45 laser periods. While this reflection occurs, 
energy fluctuations are limited by 0.04% of the total energy. 



On the other hand, when the experiment was modified by 
setting the density to 1.2nc, we observe a reflection of the 
entire electromagnetic wave by the plasma. 
We point out that these simulations have been performed using 
just the fluid part of our combined code without any PIC 
macroparticles. Still, the effect has been described correctly. 



C Target Normal Sheath Acceleration 



B. Energy conservation 

An important property we require from the new integrator is 
the conservation of the total energy of the system, comprising 
PIC macroparticles, electromagnetic fields, and the hybrid 
fluid. A very simple setup with a laser pulse being reflected 
from an overdense surface is used for this benchmark. We 
expect the total energy 

Etot = V mic\-f - 1) + ^ / (^2 ^ B^)dV 
+ / nh{^h 

JV 



l)mhC dV 



to be constant, where mi are the masses of the respective 
particle species and 7 = y/l + {pi/micY is the relativistic 
gamma factor. We denote the hybrid density by Uh and its 
gamma factor by 7^. Figure 4 shows the total energy of the 
simulation versus time, which is measured in units of laser 
periods. 



For a more realistic benchmark we model a physical setup our 
hybrid code is very suitable for: We use it for the investigation 
of the TNSA process. [24]. TNSA provides a possible way of 
laser ion acceleration out of solids by utilizing the electrostatic 
fields generated by the space charge of thermal electrons. 
The process is shown schematically in figure 5: A 10/s laser 
pulse of normalized amplitude ao = 2 is focused on a thin 
foil which can be assumed to have been pre-ionized by the 
laser. The foil consists of a bulk part of lOOOnc, a preplasma 
on its front surface, and an 80 nm thick proton layer on its 
back surface. The preplasma is modeled as a density ramp 
reaching from to 2nc over a distance of 2 laser wavelengths 
(1.6 jim) and treated entirely by the PIC method. Analogously, 
we use PIC macroparticles for the back surface protons. On the 
contrary, any attempt to describe the highly overdense main 
part of the foil as macroparticles would result in numerical 
problems. Here we use the hydrodynamic feature of H-VLPL, 
setting the hybrid density on the grid to lOOOnc. 




D. Comparison of skin depths 
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Plasma target 
Fig. 5. Schematic of the TNSA process. 
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Fig. 6. Snapshot of the TNSA benchmark simulation after 10 (left) and 380 
(right) laser periods. PIC macroparticles containing electrons are displayed 
blue, while those with protons are rendered red. One observes the thin coating 
of protons dissolving from the back of the foil in the right image. 
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Fig. 7. Spectrum of the accelerated ions in the TNSA simulation. 



The intense laser radiation creates a blow-off region in the 
front of the foil, resulting in a large cloud of hot electrons, 
which, in parts, propagates through the foil and passes the 
coating of the back surface. As the electrons leave the surface, 
a strong electrostatic field is built up, and the protons are 
pulled out of the foil and eventually accelerated to high 
energies. 

In Fig. 7, the spectrum of the accelerated ions is shown. A 
maximum energy of about 0.9 MeV is reached, which is quite 
remarkable considering the laser intensity in the setup. 
We conclude that our hybrid algorithm succeeded well and 
efficiently in treating this numerically challenging physical 
situation. 



As a further benchmark for our hybrid code we check the 
decay of a wave in an overdense plasma. According to the 
linear theory, it should scale as E{x) ~ exp(— x/(5s) in the 
plasma, where 5s = c/Juo'^ — uS^ is the skin length. Several 
simulations have been set up using a circularly polarized laser 
pulse with duration 6 A and amplitude ao = 0.01 in order to 
avoid relativistic nonlinearities. The densities of the plasma 
surfaces used for this benchmark range from X.hUc to 500nc. 
We show the decay of the wave inside the plasma for three 
densities; the agreement with the theoretical predictions up to 
densities of 500nc is very good. 



Additionally, by fitting exponentials through the measured 
field data, one can compute the skin depths of the decay. In fig- 
ure 9, the results are shown and we get an excellent agreement. 
One has to mention that even though these simulations have 
been done with a grid step of 0.05 A, the skin depths match 
remarkably well with the theory up to a density of 500nc, 
where 5^ = 0.007A. 
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Fig. 8. Snapshot of the logarithm of the fields \og{JE^ + ^1) (solid 
lines) inside the plasma for three different densities. The dashed lines show 
the theoretical prediction. 
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The setup is restricted to a 2D geometry, with the beams 
travehng perpendicularly to the x-y-plane; this is necessary 
in order to exclude two-stream instabilities. After about 3.3 
beam plasma periods 27r/a;5, with u = Y^47rn5e^/m, one 
observes a strong filamentation of the beam, and a magnetic 
field builds up. When launching the same simulation with and 
without the hybrid model, we notice that the latter succeeds 
well in describing the filamentation effect at the initial, linear 
stage. We compare the integral of the squared magnetic field 
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Fig. 9. Plot of the skin depth versus the plasma density. The blue line shows 
the grid step used in the simulation. 



E. Investigation of Weibel instabilities 

When studying the fast ignition (FI) scenario [8] in inertial 
confinement fusion, one is interested in the behaviour of the 
beam of electrons propagating into the target, particulary the 
amount of energy deposited and the shaping of the beam over 
time. Generally, electron beams running through a background 
plasma suffer from the major problem of the Weibel instability 
[25], which is a very important issue to be studied if one wants 
to understand the FI scheme. The ratio of the beam density to 
that of the background rib /rip, as well as the density gradient 
in propagation direction, is likely to influence the evolution of 
the beam, its filamentation and the increase of electromagnetic 
fields as the instability builds up. 

For low densities, roughly about lOOnc, PIC simulations can 
be carried out to perform these investigations. However, as the 
electron beam approaches the core of an ICF pellet, the density 
will exceed multiple times solid density and conventional PIC 
codes must be applied with extremely small time steps in order 
to avoid numerical instability, and thus cannot be used with 
reasonable computational effort. 

We are going to study the phenomenon of the Weibel instabil- 
ity with our new hybrid code, using standard PIC macroparti- 
cles for the electron beam and the fluid part in order to model 
the background plasma. Since H-VLPL has no restrictions for 
the hybrid densities used, we can perform such simulations 
within a moderate amount of CPU time. 
In order to obtain a proof for the physical correctness of our 
code within the linear regime, we have launched tests with 
H-VLPL comparing a classical PIC computation to a hybrid 
simulation of this setup. An electron beam with density rib 
propagates through a background plasma with rip = lOOn^. 
The momentum of the beam electrons is pb = mc with 
a thermal spread of 10 ~^mc, and the momentum of the 
background is chosen such that its current compensates for 



of the two models. At this point it has to be mentioned 
that during the nonlinear stage of the instability, the present 
version of H-VLPL will fail in describing the filamentation 
of the background plasma since it does not treat its continuity 
equation and convective term of momentum evolution. 
Additionally, the fluid plasma does not react to magnetic 
fields directly. 
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Fig. 10. Snapshot of the Weibel instability benchmark simulations with PIC 
(left) and the hybrid code (right) at 3.5, 6, and 20 beam plasma periods. We 
observe very similar behaviour, although the instability starts approximately 
0.5 period later with the hybrid model due to the lower numerical noise. 
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Fig. 11. Integral of the squared B-field over the simulation plane. During the 
linear stage of the instability, we observe an exponential growth with almost 
the same growth rate. 



Nevertheless, the behaviour, and the growth rate of the Weibel 
instabihty during the Hnear stage are accurately reproduced. 
This result indicates the applicability of H- VLPL to the Weibel 
instability scenario, and makes further investigations of the 
effect with an advanced, fully hydrodynamic hybrid code 
appear promising. 

V. Outlook 

The next step in the development of our hybrid laser plasma 
simulation system H-VLPL will be the full nonlinear hydro- 
dynamic description of the background plasma. This includes 
the continuity equation to describe the fluid transport as well 
as momentum transport equation. We are going to further 
study the physical effects mentioned above, namely the TNSA 
process and the Weibel instability, using an advanced version 
of H-VLPL, which is currently under development. 
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